#pragma once

void compute_s0_ff(colorcolorspinspin ft, complex *ff, double *p, double p2)

{       

///////////////////////////////////////////////////////////////////
// physical basis
///  S(p) = -i*pslash/p^2 sig1 + sig2/p^2 + i*gamma5/p^2*sig3   ///
// Twisted basis S^tw(p)=e(ipi/4 g5) S(p) e(ipi/4 g5)
//   S^tw(p) = -i pslash/p^2 sig1 + ig5 sig2/p^2 - 1/p^2*sig3  
///////////////////////////////////////////////////////////////////
        int idir;
        double pmu;

        ff[0]=ff[1]=ff[2]=0.;

        /////////////
        // Sigma1 //
        //////////////
       for (int mu=1; mu<5; mu++){
             if (mu==4) idir=0;
             else idir=mu;
             pmu=p[idir];
             for (int ic=0; ic<3; ic++)
                for (int is=0; is<4; is++)
	        {
		ff[0]=ff[0]+I*(base_gamma[mu]).entr[is]*ft[ic][ic][(base_gamma[mu]).pos[is]][is]*pmu/12.;
                }
       }

       /////////////
       // Sigma2 //
       ////////////
       for (int ic=0; ic<3; ic++)
          for (int is=0; is<4; is++)
          {
		ff[1]=ff[1]+ft[ic][ic][is][is]/12.;
          }
        ff[1]=ff[1]*p2;

        /////////////
        // Sigma3 //
        ////////////
        for (int ic=0; ic<3; ic++)
          for (int is=0; is<4; is++)
          {
	   ff[2]=ff[2]-I*(base_gamma[5]).entr[is]*ft[ic][ic][(base_gamma[5]).pos[is]][is]/12.;
          }
          ff[2]=ff[2]*p2;
}

void compute_s0tw_ff(colorcolorspinspin ft, complex *ff, double *p, double p2, double nn, int r)

{

///////////////////////////////////////////////////////////////////
// physical basis
///  S(p) = -i*pslash/p^2 sig1 + sig2/p^2 + i*gamma5/p^2*sig3   ///
// Twisted basis S_{+}^tw(p)=e(ipi/4 g5) S(p) e(ipi/4 g5)
//               S_{-}^tw(p)=e(-ipi/4 g5) S(p) e(-ipi/4 g5)
//   S_{+}^tw(p) = -i pslash/p^2 sig1 + ig5 sig2/p^2 - 1/p^2*sig3  
//   S_{-}^tw(p) = -i pslash/p^2 sig1 - ig5 sig2/p^2 + 1/p^2*sig3  
///////////////////////////////////////////////////////////////////

	int pm;
	if (r==0) pm=1;
	else if (r==1) pm=-1;

        int idir;
        double pmu;

        ff[0]=ff[1]=ff[2]=0.;

                     /////////////
                     // Sigma1 //
                     //////////////
                     for (int mu=1; mu<5; mu++){

                                if (mu==4) idir=0;
                                else idir=mu;

                                pmu=p[idir];
                                for (int ic=0; ic<3; ic++)
                                  for (int is=0; is<4; is++)
                                   {
					ff[0]=ff[0]+I*(base_gamma[mu]).entr[is]*ft[ic][ic][(base_gamma[mu]).pos[is]][is]*pmu/12.;

                                   }
                       }

                   /////////////
                   // Sigma2 //
                   ////////////
                  for (int ic=0; ic<3; ic++)
                     for (int is=0; is<4; is++)
                     {
                                ff[1]=ff[1]-pm*I*(base_gamma[5]).entr[is]*ft[ic][ic][(base_gamma[5]).pos[is]][is]/12.;
                     }
                     ff[1]=ff[1]*p2;

                  /////////////
                  // Sigma3 //
                  ////////////
                 for (int ic=0; ic<3; ic++)
                   for (int is=0; is<4; is++)
                   {
                           ff[2]=ff[2]-pm*ft[ic][ic][is][is]/12.;

                   }
                 ff[2]=ff[2]*p2;
}
void compute_s0_inv_ff(colorcolorspinspin ft, complex *ff,double *p, double p2, double nn)

{

///////////////////////////////////////////////////////////////////
///  S(p) = -i*pslash sig1 + sig2 - i*gamma5*sig3   ///
///////////////////////////////////////////////////////////////////


        int idir;
        double pmu;


        ff[0]=ff[1]=ff[2]=0.;

                     /////////////
                     // Sigma1 //
                     //////////////
        for (int ic=0; ic<3; ic++)
             for (int is=0; is<4; is++){
                     for (int mu=1; mu<5; mu++){

                                if (mu==4) idir=0;
                                else idir=mu;
                                pmu=p[idir];

                                if (fabs(pmu)>1e-10){
                                ff[0]=ff[0]-I*(base_gamma[mu]).entr[is]*ft[ic][ic][(base_gamma[mu]).pos[is]][is]/12./pmu/nn;
                                }

              }
         }

                   /////////////
                   // Sigma2 //
                   ////////////
                  for (int ic=0; ic<3; ic++)
                     for (int is=0; is<4; is++)
                     {
                                ff[1]=ff[1]+ft[ic][ic][is][is]/12.;
                     }

                  /////////////
                  // Sigma3 //
                  ////////////
                 for (int ic=0; ic<3; ic++)
                   for (int is=0; is<4; is++)
                   {
                           ff[2]=ff[2]-I*(base_gamma[5]).entr[is]*ft[ic][ic][(base_gamma[5]).pos[is]][is]/12.;
                   }
}
void compute_s0tw_inv_ff(colorcolorspinspin ft, complex *ff,double *p, double p2, double nn, int r)

{

///////////////////////////////////////////////////////////////////
///  S_{+-}tw(p) = -i*pslash sig1 +- ig5 sig2 +- sig3   ///
///////////////////////////////////////////////////////////////////


        int pm;
        if (r==0) pm=1;
        else if (r==1) pm=-1;


        int idir;
        double pmu;


        ff[0]=ff[1]=ff[2]=0.;

                     /////////////
                     // Sigma1 //
                     //////////////
        for (int ic=0; ic<3; ic++)
             for (int is=0; is<4; is++){
                     for (int mu=1; mu<5; mu++){

                                if (mu==4) idir=0;
                                else idir=mu;
                                pmu=p[idir];
				if (fabs(pmu)>1e-10){
                                ff[0]=ff[0]-I*(base_gamma[mu]).entr[is]*ft[ic][ic][(base_gamma[mu]).pos[is]][is]/12./pmu/nn;
				}

              }
         }

                   /////////////
                   // Sigma2 //
                   ////////////
                  for (int ic=0; ic<3; ic++)
                     for (int is=0; is<4; is++)
                     {
	                           ff[1]=ff[1]-pm*I*(base_gamma[5]).entr[is]*ft[ic][ic][(base_gamma[5]).pos[is]][is]/12.;

                     }

                  /////////////
                  // Sigma3 //
                  ////////////
                 for (int ic=0; ic<3; ic++)
                   for (int is=0; is<4; is++)
                   {
                              ff[2]=ff[1]+pm*ft[ic][ic][is][is]/12.;

                   }
}

